################################################################################
# Created By:Pietryka
# Creation Date:  2016-08-22
# Purpose: Fit models and format Table 2
# Questions: mpietryka@fsu.edu
################################################################################

# PREAMBLE  ====================================


# LOAD PACKAGES --------------
library(tidyverse)


# LOAD DATA --------------

# data from 'CPS-1B-Clean-Participation.R'
load("Data/CPS-1B-Clean-Participation.Rdata")


# CLEAN DATA   ====================================

cps_combined <- cps_sub_10  %>%
  tbl_df()  %>%
  mutate(
    # EMPLOYMENT
    work_employed =  if_else(PEMLR == 1, 1, 0, missing = NA_real_),
    work_retired  =  if_else(PEMLR == 5, 1, 0, missing = NA_real_),
    # RELATIONSHIP STATUS
    marry_married  =  if_else(PEMARITL == 1, 1, 0, missing = NA_real_),
    marry_sepordiv  =  if_else(PEMARITL %in% c(2, 4, 5) ,
                               1,
                               0,
                               missing = NA_real_),
    marry_widow  =  if_else(PEMARITL == 3, 1, 0, missing = NA_real_),
    # EDUCATION
    educ_hs = if_else(educ == 39, 1, 0, missing = NA_real_),
    educ_somecol = if_else(educ > 39 & educ < 43,
                           1,
                           0,
                           missing = NA_real_),
    educ_colgrad = if_else(educ == 43, 1, 0, missing = NA_real_),
    educ_advanced = if_else(educ > 43 & educ <= 46,
                            1,
                            0,
                            missing = NA_real_),
    # RACE
    race_black = if_else(PTDTRACE == 2, 1, 0, missing = NA_real_),
    race_asian = if_else(PTDTRACE == 4, 1, 0, missing = NA_real_),
    race_other = if_else(PTDTRACE %in% c(3, 5, 6:21), 1, 0, missing = NA_real_),
    # INCOME
    famincquart = ntile(famincome, 4)
  )



# PREPARE MODELS   ====================================


# EXPLANATORY VARIABLE AND CONTROL NAMES
righthand_names <- c(
  "prekkids",
  "age",
  "female",
  "student",
  "work_employed",
  "work_retired",
  "educ_hs",
  "educ_somecol",
  "educ_colgrad",
  "educ_advanced",
  "marry_married",
  "marry_sepordiv",
  "marry_widow",
  "factor(famincquart)",
  "race_black",
  "race_asian",
  "race_other"
)


outcome_names[which(outcome_names == "discuss_pol")] <- "as.factor(discuss_pol)"
formulae <- lapply(outcome_names, function(x){
  paste0(x, " ~ ", paste(righthand_names , collapse = " + "))  %>%
    as.formula()
})


form_mod <- data_frame(
  formula = formulae,
  mod_type = c("ols", rep("logit", length(outcome_names) - 1))
)

fit_model <- function(the_formula, mod_type) {
  if (mod_type == "ols") {
    MASS::polr(the_formula, data =  cps_combined, Hess = TRUE)
  } else{
    if (mod_type == "logit") {
      glm(the_formula, data =  cps_combined, family = "binomial")
    }
  }
}

# FIT MODELS   ====================================
library(texreg)

cps_mods <- form_mod  %>%
  mutate(mod = map2(formula, mod_type, fit_model),
         outcome_var = outcome_names,
         outcome_lab = c("1. Discuss", "2. Contact", "3. Community", "4. Civic", "5. Officer"),
         mod_extract = map(mod, extract,  include.thresholds = TRUE)
  )




var_map <- list(
  prekkids = "Parent of Young Child (0 = No; 1 = Yes)",
  age = "Age (in Years)",
  female = "Female (0 = No; 1 = Yes)",
  student = "Student (0 = No; 1 = Yes)",
  work_employed = "Employed (0 = No; 1 = Yes)",
  work_retired = "Retired (0 = No; 1 = Yes)",
  educ_hs = "HS Graduate (0 = No; 1 = Yes)",
  educ_somecol = "Some College (0 = No; 1 = Yes)",
  educ_colgrad = "College Graduate (0 = No; 1 = Yes)",
  educ_advanced = "Advanced Degree (0 = No; 1 = Yes)",
  marry_married = "Married (0 = No; 1 = Yes)",
  marry_sepordiv = "Separated/Divorced (0 = No; 1 = Yes)",
  marry_widow = "Widowed (0 = No; 1 = Yes)",
  `factor(famincquart)2` = "In Second Quartile (0 = No; 1 = Yes)",
  `factor(famincquart)3` = "In Third Quartile (0 = No; 1 = Yes)",
  `factor(famincquart)4` = "In Fourth Quartile (0 = No; 1 = Yes)",
  race_black = "Black (0 = No; 1 = Yes)",
  race_asian = "Asian (0 = No; 1 = Yes)",
  race_other = "Other (0 = No; 1 = Yes)",
 `(Intercept)` = "0|1 (Intercept)",
 `0|1`       = "0|1 (Intercept)",
 `1|2`       = "1|2",
 `2|3`       = "2|3",
 `3|4`       = "3|4"
)




group_list <- list(
  " " = 1:4,
  "Employment Status (Reference = Unemployed/Not in Labor Force)" = 5:6,
  "Education (Reference = Less than HS Graduate)" = 7:10,
  "Relationship (Reference = Never Married)" = 11:13,
  "Income (Reference = In First Quartile)" = 14:16,
  "Race (Reference = White)" = 17:19,
  "Thresholds" = 20:23
)

screenreg(cps_mods$mod_extract,
          custom.coef.map = var_map,
          groups = group_list)


screenreg(cps_mods$mod_extract)

# ------
htmlreg(cps_mods$mod_extract ,
        caption = "Table 2. Regressions of Non-voting Participation on Parenthood Status and Controls using CPS Data.",
        caption.above = TRUE,
        stars = 0.05,
        groups = group_list,
        custom.coef.map = var_map,
        custom.model.names = cps_mods$outcome_lab,
        single.row = TRUE,
        custom.note = "* p < 0.05. Source: CPS 2010 Civic Engagement Supplement",
        file = "Tables/table-CPS-participation.html"
)

